Nonlinear wave-packet dynamics in a disordered medium 
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In this article we develop an effective theory of pulse propagation in a nonlinear and disordered 
medium. The theory is formulated in terms of a nonlinear diffusion equation. Despite its apparent 
simplicity this equation describes novel phenomena which we refer to as "locked explosion" and 
"diffusive" collapse. In this sense the equation can serve as a paradigmatic model, that can be 
applied to such distinct physical systems as laser beams propagating in disordered photonic crystals 
or Bose- Einstein condensates expanding in a disordered environment. 



In recent years, novel experimental techniques made 
possible first observations of wave-packets evolving in the 
presence of random scatterers and nonlinearities. In a 
number of optical experiments, a laser beam was sent 
into a nonlinear optical medium with a random refrac- 
tive index, and the beam profile in the transverse di- 
rection (s) was monitored on the opposite side of the 
sample [3, Q. In a second class of experiments, atoms 
forming a Bose-Einstein condensate were released from 
a trap and subjected to a disorder potential during the 
expansion 0, 0, @, 0, 0]- The experiments were in- 
spired by the idea that in these setups, unlike for trans- 
port experiments in electronic systems, one can visual- 
ize the phenomenon of Anderson localization, whereby 
a wave-packet or quantum particle is confined within 
a finite volume as a result of multiple scattering on a 
random potential. The evolution of the injected wave- 
packet in both experiments can be described by the 
non-linear Schrodinger equation (NLSE), in the context 
of atomic Bose-Einstein condensates referred to as the 
Gross-Pitaevskii Equation (GPE). This equation differs 
from the linear Schrodinger equation by an additional 
cubic term and is used as an exemplary description for 
nonlinear waves. The nonlinear ity is a consequence of 
interactions between particles in the case of atomic con- 
densates and of a change in the refractive index in re- 
sponse to the electric field (Kerr effect) in the case of 
laser beams. Motivated by these experiments we derive, 
starting from the GPE/NLSE, a kinetic equation that 
describes the evolution of an injected wave-packet in a 
weakly disordered nonlinear medium in two dimensions. 
Analysis of this equation reveals a rather nontrivial pic- 
ture: Irrespective of the sign of the nonlinearity the mean 
square radius of the wave-packet changes linearly in time, 
dt (r 2 ) oc E to t, where E tot is the total energy of the wave- 
packet. For a repulsive nonlinearity the initial change 
of the profile displays features of an explosion, although 
the overall size of the wave-packet is growing slowly as 
in ordinary diffusion. For an attractive nonlinearity, the 
radius can either grow or decrease, depending on the sign 
of E tot . In particular, for E tot < we predict a slow "dif- 
fusive" collapse as the radius of the wave-packet shrinks 
towards zero. 



In this paper we will mostly use the language related 
to the GPE, but also indicate below how to translate to 
a language more suitable for optical experiments. The 
GPE describes the evolution of a macroscopic wave- 
function ^ Q: 



1 



(1) 



= - ^ V 2 ^(r, t) + u(r)tf (r, t) + A|tt(r, t)| 2 *(r, f), 

where we set ft = 1. For positive (negative) A this equa- 
tion contains a repulsive (attractive) self-consistent po- 
tential A|\£(r, t)\ 2 . This corresponds to a nonlinearity 
of the de-focusing (self-focusing) type. The static dis- 
order potential u(r) is the source of randomness in the 
equation. For simplicity we choose for our calculation a 
Gaussian white noise potential with correlation function 
(u(r)u(r f )) = S(r — r / )/(mr) [for a discussion of averag- 
ing for speckle potentials see, e.g., Ref. [lOj]]. The angu- 
lar brackets denote averaging over disorder configurations 
and r is the scattering time. 

The NLSE used in optics is derived in the so-called 
paraxial approximation , and thus describes the evo- 
lution of the smooth envelope of the electric field. The 
main propagation direction of the laser beam, say the z- 
direction, plays the role of time in the NLSE. In this 
sense, the disorder potential which results from ran- 
dom variations of the refractive index is static when 
it is z-independent. For example, the two-dimensional 
transverse evolution of a pulse is studied in a three- 
dimensional sample. The intensity of the beam is pro- 
portional to |^(r, z)\ 2 . In the NLSE, the mass m in the 
GPE is replaced by the wave vector k = lj/c, where uo 
is the frequency of the carrier wave and c the velocity of 
light in the medium. 

For a condensate released from a confining harmonic 
oscillator potential, as is typical for experiments on cold 
atomic gases, the GPE without disorder can be solved 
exactly [12|, ll3[ . During an initial stage the potential en- 
ergy originating from the nonlinearity is almost entirely 
converted into kinetic energy. This period of violent ac- 
celeration is followed by a second stage, during which 
the nonlinearity is no longer essential. Expansion in the 
presence of disorder in the two-dimensional case was re- 
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cently addressed in reference [14]. In this paper it has 
been assumed that for the repulsive nonlinearity an ini- 
tial ballistic stage is not affected by disorder, while the 
subsequent diffusive expansion is not affected by the non- 
linearity, thereby separating the two effects. In contrast, 
we are interested here in the interplay of disorder and 
nonlinearity, both attractive and repulsive. This is espe- 
cially interesting in two dimensions, as it is known that 
for linear wave propagation and weak disorder there is an 
extended diffusive regime preceding localization on larger 
length scales. 



n(r,£) 



It turns out that a perturbation theory organized in 
powers of the interaction constant A is not well suited 
for the nonlinear problem, and hence a non-perturbative 
approach is required. Since the system is far out of equi- 
librium, we choose to work with a kinetic equation. The 
derivation of the kinetic equation proceeds as follows. We 
use methods of classical statistical field theory to derive 
a functional integral expression for the disorder averaged 
density 15, [l6j]. The formalism involves a doubling of 
the degrees of freedom, similar to the Keldysh or closed- 
time-path approaches for quantum systems [l6j], where 
two fields are introduced on forward and backward time- 
contours. Instead of averaging over a statistical ensemble 
in the initial state, we assume that the wave-function at 
the initial time is known. Averaging is performed over 
disorder configurations. Scattering on impurities is in- 
cluded on the level of the self-consistent Born approx- 
imation. While interference (weak localization) correc- 
tions are not covered by this approximation, it allows 
for a consistent description of diffusion in the presence of 
nonlinearity we are focusing on in this paper. The nonlin- 
earity is treated by introducing a self-consistent potential 
$(r, t). In this way it is possible to include interaction 
effects in a non-perturbative way, which is crucial for the 
problem at hand. To obtain the kinetic equation for the 
density in the diffusive limit, we assume that the initial 
wave-function sets a momentum scale po characterizing 
the main part of the momentum distribution, so that 
the weak disorder condition pol ^> 1 is fulfilled, where 
I — p^Tjm is the mean free path. We further assume 
that the density varies smoothly on scales of /, in partic- 
ular that the size of the condensate is much larger than 
the mean free path. Both of these conditions can be met 
simultaneously. The phase of which is related to the 
momentum, may change rapidly, while the amplitude, 
which determines the density, may vary smoothly. Even 
if the density does not satisfy the smoothness condition 
initially, it is natural to expect that in the case of an 
expansion it will become sufficiently smooth after some 
time. The derivation of the kinetic equation will be pre- 
sented elsewhere [17j . 

Starting from Eq. (pQ), the outlined steps lead to the 
following equation for the density evolution in the diffu- 



FIG. 1: Diagrammatic illustration. The injection process 
takes place on the left hand side and the time arrow points 
to the right. The solid lines are disorder averaged retarded 
(right-arrow) and advanced (left-arrow) Green's functions. 
The dashed lines describe the result of averaging over dis- 
order. They form ladder diagrams, which graphically repre- 
sent a diffusion process. The particular way of averaging over 
disorder is justified in the limit SkinT ^> 1, where Skin is de- 
fined below Eq. J2}. The wavy lines account for the effective 
interaction induced by the nonlinearity in Eq. JT]). 



sive regime 

d t n(r, t, e) - V(ZVtf Vn(r, t, e)) + d t tf (r, t)d e n(r, t, e) 
= 5(t)F(e-#(r,0),r) (2) 

where F(e,r) = J[d 2 qd 2 p/ (2tt) 4 } F(p,q) exp(iqr) 2n5(s- 
e p ), and F(p,q) = ^ (p + q/2)^(p - q/2) is deter- 
mined by the initial wave function \I>o; £p = p 2 /(2m) 
is the kinetic energy, and D £ = er/m is the diffusion 
coefficient. The equation should be supplemented 
with the self-consistency relation for the potential 
i?(r, t) = 2An(r, t), where n(r, t) = j ds/(27r) h(r,t,e). 
Despite its apparent simplicity it is a rather complicated 
nonlinear integro-differential equation. The kinetic 
equation effectively sums an infinite series of diagrams 
of the type shown in Fig. [TJ It is a peculiarity of 
the perturbation theory for a classical field equation 
such as the GPE that no closed loops arise, making it 
quite distinct from the related problem in interacting 
electron systems. The relation between certain blocks 
appearing in the diagrammatic perturbation theory 
and the corresponding terms in the kinetic equation is 
visualized in Fig. [2j 

The physics described by Eq. ([2]) is essentially classi- 
cal. Imagine first that the potential $ does not depend 
on time. Consider now a particle diffusing with total en- 
ergy £ on the background of a smoothly varying potential 

If scattering events are frequent enough, the diffusion 
coefficient is determined by the kinetic energy e p = e — d 
that varies locally in space. If the potential additionally 
varies in time, the particle may change its total energy. 
If on the other hand the potential depends on time only, 
the kinetic energy does not change. Therefore, it is ex- 
pected that a purely time-dependent potential has no ef- 
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FIG. 2: Relation between the diagrammatic representation 
and the kinetic equation. The block shown in (a) gives rise to 
the term V(D e -#Vn(r, t, e)) and the one in (b) to the term 
dt$(r, t)d £ n(r, £, e) in the kinetic equation, Eq. J2]). 



feet on the density. This observation is related to the fact 
that in the original GPE a purely time dependent poten- 
tial V(t) may be removed by a suitably chosen gauge- 
transformation, *(r,t) -> #(r,t)exp(-i/* o dt'y (£')), 
that does not affect the density |^(r, t)\ 2 . Indeed, we can 
make this point obvious in Eq. ([2]) by shifting the energy 
variable so that it will correspond to the kinetic energy 
instead of the total energy, n(r, e, t) = n(r, e + #(r, t), t). 
Expressed in the new coordinates the equation reads 



d t n(r,e,i) 

- V P - V r tf(r,£)<9 £ 



V r - V r tf(r,t)<9 £ n(r,e,t) 
(3) 



An equation for the density n(r, £) can be obtained by 
integrating Eq. (|3j) in e 

d t n(r,t) - -V 2 (e(r,t) + An 2 (r,£)) = n(r,0), (4) 
m 

where e(r,i) = J ds/(27r) en(r,t,e) = £&m(r, £)n(r, £). 
It can be written in a more compact form, dtn — 
V 2 (D e ffn) = 5(t)n, when defining an effective space 
and time-dependent diffusion coefficient D e ff = (skin + 
Xn)r/m. The apparent simplicity of these equations is 
however deceiving. They are not closed equations for the 
density evolution, since the kinetic energy s generally de- 
pends on the strength of the nonlinearity and needs to 
be determined separately via Eq. ([3|). Nevertheless, we 
arrive in this way at the conceptually important result 
that in the diffusive regime the nonlinearity effectively 
introduces a density dependence into the diffusion coef- 
ficient. 

It seems clear that a closed form solution of the non- 
linear equations for arbitrary initial conditions cannot be 
found. In order to make progress we will rely on two ap- 
proaches: the use of conservation laws and the study of 
solvable limiting cases. When combined, they will enable 
us to arrive at a qualitative picture both for repulsive and 
attractive nonlinearity. 



First we briefly discuss the linear case, $ = 0. In the 
absence of nonlinearity, n(r, e, t) evolves independently 
for each energy e. In this limit, Eq. (|3]) has the obvious 



solution n(r, e, t) 



e(t) 



4tt£> 



i e" 



-(r- ri )7(4D £ £) 



F{e,n). 



Here, for each energy e diffusion is determined by the 
corresponding diffusion coefficient D e , and should be 
weighted according to the energy distribution in the 
injected wave-packet [18] . This case was discussed in 
Ref. [3. 

Next we turn to the nonlinear case. We will make use 
of the conservation laws for particle number and energy. 
By integrating Eq. (j4]) over r, we immediately obtain that 
the particle number (or normalization) J dr n(r, t) = N 
is fixed in time. An equation for e(r, t) can be derived 
by first multiplying Eq. ([3j) by e before integrating in 
this variable. Then by combining the equation for n(r, t) 
with the equation for e(r, t) we find that the en ergy 
Etot = J dr (e(r,i) + An 2 (r, t)) is constant in time [19j . 
Surprisingly, this conservation law completely determines 
the time evolution of the mean radius squared of the 
wave-packet, (r 2 ) = J dr r 2 n(r,t)/N. Indeed, multi- 
plying Eq. (|4j) by r 2 and subsequently integrating in r 
one obtains that dt (r 2 ) = 4D £tot , where e tot = E tot /N. 
The linear dependence of the mean square radius on time 
during the whole evolution is guarded by energy conser- 
vation. This is one of the central results of this paper. 
When compared to the linear case, the effective diffusion 
coefficient D £tot is reduced for attractive and enhanced 
for repulsive nonlinear it ies. 

In the following we discuss more specifically the re- 
pulsive and attractive cases. For the repulsive nonlinear 
case it is instructive to consider a situation in which the 
second term on the RHS of Eq. (|4|) dominates. The equa- 
tion dtn = V 2 n 2 , which one obtains after simple rescal- 
ing, is an example of the famous porous medium equa- 
tion (PME) 0. For the 2d case the solution describing 
the evolution of a delta- function pulse MS(r) is given by 
n(r,t) = (C-r 2 /(16t 1 / 2 ))/t 1 /2 ? where C 2 = M/(8tt) 



221 ]. This solution is often referred to as Barenblatt's 
solution. It conserves the normalization J dr n(r, t) = M 
but, unlike ordinary diffusion, it is nonzero only in a fi- 
nite region of space. The special importance of Baren- 
blatt's solution in the theory of the PME is related to the 
fact that, roughly speaking, any solution starting from a 
sufficiently benign initial pulse with weight M is eventu- 
ally well- approximated by Barenblatt's solution with the 
same weight, a statement known as the nonlinear central 
limit theorem for the PME 0. 

For Barenblatt's solution, the mean radius squared 
evolves as (r 2 ) ex t 1 / 2 and the density at r = drops 
as n(0, t) oc t~ x l 2 . At short times this solution describes 
a much faster " explosive" evolution than the source- type 
solution of the diffusion equation, for which (r 2 ) oct and 
n(0,t) ex t -1 . At first sight there seems to be a contra- 
diction. If one injects a bell-shaped pulse with a large 



potential energy, it appears that the potential part of 
the effective diffusion coefficient D e ff = (skin + Xn)r/m 
dominates. Therefore, naively, one would assume that 
the initial evolution is " explosive" , while our exact result 



-4D £tot t rules out this possibility. This puzzle 



can be resolved in the following way. The explosion takes 
place only in the central part of the density distribution, 
which has only a small weight when calculating (r 2 ). 
Right in the center, for r = 0, Eq. (|4]) can be written 
as d t n — ^ [V 2 £ + 2AnV 2 n] for t > 0; we consider here 
a rotationally symmetric distribution with Vn(0, t) = 
and V 2 n(0,0) < 0. For sufficiently large An, the po- 
tential part is dominant and leads to a fast initial de- 
crease of the density before either AnV 2 n becomes small 
or V 2 £ becomes positive as a consequence of the outward- 
flow of the kinetic energy. Away from the center, where 
the density and correspondingly the term 2\nV 2 n are 
small, Eq. (@|) takes the form d t n w ^ [V 2 £ + 2A(Vn) 2 ] 
for t > 0. For the PME it is the second term that de- 
termines the propagation of the boundary. For Eq. (|4]), 
however, the large kinetic energy outside the center leads 
to V 2 £ < for intermediate distances, and this prevents 
the term 2A(Vn) 2 from dominating. It is therefore an 
inversion of the distribution of kinetic energy compared 
to that of the density that does not allow for an explo- 
sive expansion and leads to a linear dependence of (r 2 ) 
on time. A sketch of a typical density evolution expected 
for this kind of " locked explosion" is presented in the first 
line of Fig. 3. 

We now discuss general features of wave-packet dy- 
namics in the disordered and nonlinear medium. For an 
expanding wave-packet, i.e. E tot > 0, the overall poten- 
tial energy related to the nonlinearity is converted into 
kinetic energy. As a result, the total kinetic energy in- 
creases in the repulsive case and decreases in the attrac- 
tive case. Correspondingly, during the course of the ex- 
pansion localization effects can be expected to be weak- 
ened for repulsive nonlinearity and enhanced for attrac- 
tive nonlinearity. In particular, for an attractive (self- 
focusing) nonlinearity the slowing down and eventual lo- 
calization of the injected pulse (not considered here) oc- 
curs at smaller distances than in the linear case as ob- 
served in the experiment [l], Q . 

The attractive case is richer than the repulsive one (see 
Fig. 3), because the total energy may also be negative, 
Etot < 0- Then the mean radius squared would become 
equal to zero after a finite time. This corresponds to a 
celebrated phenomenon in nonlinear physics, the collapse 
[24L [25t . Here it is realized for the diffusive system. To 
the best to our knowledge, this "diffusive" collapse has 
not been discussed in the literature. Since our reasoning 
is based on a diffusive kinetic equation and thus assumes 
frequent scattering, the linear decrease of (r 2 ) only holds 
as long as the radius of the cloud exceeds the mean free 
path [261 ]. 

Even for E tot > the collapse can play a role when 
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FIG. 3: Time-evolution of the disorder averaged den- 
sity. The time- evolution of the disorder averaged density 
n is sketched for the three relevant cases, starting from the 
same initial distribution at t = 0. According to the rela- 
tion dt (r 2 ) = Etotr/m, the radius grows for Etot > and 
shrinks for Etot < 0. The condition Etot > can be realized 
for positive (first line) or negative potential $ = An (second 
line), while Etot < can be realized only for negative $ (third 
line) . In the case of repulsion a rapid drop of the density dis- 
tribution is expected in the center for large i9 = An. This case 
is referred to as "locked explosion" in the text, since despite a 
rapid change of the density profile, (r 2 ) grows only linearly in 
t. For Etot < a "diffusive" collapse is expected. In the cases 
with $ < a fragmentation of the cloud may occur depending 
on the initial conditions. In general, it is important that the 
evolution is not only determined by the initial density, but 
also by the initial energy distribution. 



the nonlinearity is attractive, if part of the cloud has a 
negative energy, while the remaining part expands. As 
a result one can expect a fragmentation of the cloud. If 
a part of the cloud with a positive but small energy lags 
behind, this fragment may have a strong tendency to lo- 
calize. One may expect that this kind of localized or col- 
lapsing fragment generically remains from an expanding 
cloud with E tot > but attractive nonlinearity. 

To conclude, we found that the the nonlinear diffu- 
sion equation discussed in this paper contains rich physics 
that invites further numerical, analytical and experimen- 
tal investigations. 
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